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ABSTRACT 

Observational studies show that the probability of finding gas giant planets 
around a star increases with the star's metallicity. Our latest simulations of 
disks undergoing gravitational instabilities (GI's) with realistic radiative cooling 
indicate that protoplanetary disks with lower metallicity generally cool faster 
and thus show stronger overall Gl-activity. More importantly, the global cooling 
times in our simulations are too long for disk fragmentation to occur, and the 



- 2 - 



disks do not fragment into dense protoplanetary clumps. Our results suggest 
that direct gas giant planet formation via disk instabilities is unlikely to be the 
mechanism that produced most observed planets. Nevertheless, GI's may still 
play an important role in a hybrid scenario, compatible with the observed metal- 
licity trend, where structure created by GI's accelerates planet formation by core 
accretion. 

Subject headings: accretion, accretion disks — hydrodynamics — instabilities — 
planetary systems: formation — planetary systems: protoplanetary disks 

1. INTRODUCTION 

The past decade has seen the discovery of over 150 exoplanets (http:/ /www. obspm.fr/planets). 
One statistical trend that has emerged from these data is that the probability of finding a 
gas giant planet around a star with current techniques increases with the host star's metal- 
hcity (Santos et al. 2001; Fischer & Valenti 2005). As shown by Fischer & Valenti (2005), 
the high metal content of planet host stars seems to be primordial. Therefore, this trend, if 
real (Sozzetti et al. 2005), indicates that short-period gas giant planets are more likely to 
occur in metal-rich than in metal-poor protoplanetary disks. The two contending theories 
for gas giant formation are core accretion plus gas capture (Pollack et al. 1996) and disk 
instability (Boss 2002; Mayer ct al. 2004). Calculations show that the metallicity relation 
can be explained within the framework of core accretion (e.g., Ida & Lin 2004; Kornet et al. 
2005). For disk instability. Boss (2002) finds that, in his three-dimensional hydrodynamics 
disk simulations with radiative cooling, clump formation by disk instability occurs for all 
metaUicities over the range 0.1 to 10 Zq, due to rapid coohng by convection (Boss 2004), 
and he attributes the abundance of short period gas giants around high metallicity stars to 
migration (Boss 2005), a mechanism also invoked by Ida & Lin (2004) to explain part of 
the metallicity correlation. By contrast, Mejfa (2004), who uses a somewhat more sophisti- 
cated treatment of radiative boundary conditions, finds much longer cooling times and no 
fragmentation into dense clumps in her disk instability simulations. Here we report results 
of new disk calculations based on Mejia's methods in which the opacity is varied by using 
different metaUicities and grain sizes. Even over a much narrower range of metaUicities than 
considered by Boss (2002), we find that the strength of the GI's does vary noticeably and 
that disk fragmentation is not seen for any metallicity or grain size tested. 
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2. 



METHODS 



We conduct protoplanetary disk simulations using the Indiana University Hydrodynam- 
ics Group code (see Pickett et al. 1998, 2000, 2003; Mejfa 2004; Mejfa et al. 2005), which 
solves the equations of hydrodynamics in conservative form on a cylindrical grid (r, 0, z) to 
second order in space and time using an ideal gas equation of state. Self-gravity and shock 
mediation by artificial bulk viscosity are included. Reflection symmetry is assumed about 
the equatorial plane, and free outflow boundaries are used along the top, outer, and inner 
edges of the grid. 

We adopt the treatment of radiative physics detailed in Mejfa (2004) with few modifl- 
cations. Let tr be the optical depth, deflned by using the Rosseland mean opacity measured 
vertically down from above. Energy flow in cells with tr > 2/3 is calculated in all three 
directions using flux-limited diffusion (Bodenheimer et al. 1990). Cells with < 2/3, in the 
disk atmosphere and in the outer disk, cool radiativcly using an optically thin LTE emis- 
sivity. Atmosphere heating by high-altitude shocks and upward moving photons from the 
photosphere are included. In this paper, we also assume that an external envelope heated 
by the star (Natta 1993; D'Alessio, Calvet & Hartmann 1997) shines vertically down on 
the disk. This IR irradiation is characterized by a black body flux with a temperature Tirr- 
The optically thick and thin regions are coupled, over one or two cells, by an Eddington-like 
grey atmosphere flt that deflnes the boundary flux for the diffusion approximation. The 
opacities and molecular weights for a solar composition are from D'Alessio et al. (2001), 
with a power-law grain size distribution of n(a) ~ a~^'^ ranging from 0.005 //m to a largest 
grain size a^ax thai can be varied. To model variations in metallicity Z, the mean opacities 
are multiplied by a factor = Z/Zq, as was done by Boss (2002). Tests of our radiative 
scheme for a vertically stratified gas layer with a constant gravity, a constant input flux at 
the base, and a grey opacity show relaxation to an Eddington-like solution with the correct 
flux from the photospheric layers. 



The initial axisymmetric model for all the calculations is the same as that used by 
Mejfa et al. (2005). The central star is 0.5 Mq, and the nearly Keplerian disk of 0.07 Mq 

has a surface density E(r) oc r~°'^ from 2.3 AU to 40 AU. The initial grid has (256, 128, 
32) cells in (r, 0, z) above the midplane. When the disk expands at the onset of GI's, the 
grid is extended radially and vertically. The initial minimum value of the Toomre stability 



3. 



SIMULATIONS 



3.1. 



Initial Model and the Set of Simulations 
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parameter Q is about 1.5, and so the disk is marginally unstable to GI's. The initial model 
is seeded with low-amplitude random density noise. We use Ti^j. = 15 K, which is lower 
than the 50 K assumed in Boss (2002) because our larger and less massive disk is mostly 
stabilized by Tj„. = 50 K. In this paper, we present simulations with four metallicities Z = 
1/4 Zq (one-quarter solar metallicity), 1/2 Zq, Zq, and 2 Zq. The 1/4 Zq simulation was 
started from the 1/2 Zq disk after 13.0 outer rotation periods (ORPs) of evolution, to save 
computing resources. Here 1 ORP (about 250 yrs) is the initial rotation period at 33 AU. 
The varied metallicity cases use a maximum grain size a^ax — 1a*™ in the dust opacities. 
An additional simulation with Umax — 1 mm and Z = Z© is conducted to explore the effects 
of grain growth. 



3.2. Results 

The current calculations resemble those presented in Mejia (2004) and Meji'a et al. 
(2005). The disks remain fairly axisymmetric until a burst phase of rapid growth in non- 
axisymmetric structure. Subsequently, the disks gradually transition into a quasi-steady 
asymptotic phase, where heating and cooling are in rough balance, and average quantities 
change slowly (see also Lodato & Rice 2005). Table 1 summarizes some of the results. In the 
tabic, Duration refers to the simulation length measured in ORPs, ti is the time in ORPs 
at which the burst phase begins, t2 is the approximate time in ORPs when the simulation 
enters the asymptotic state, {A) is a time-averaged integrated Fourier amplitude for all non- 
axisymmetric structure (see below) , tcooi is the final global cooling time obtained by dividing 
the final total internal energy of the disk by the final total luminosity, and Thin% is the 
percentage of disk volume that is optically thin during the asymptotic phase. 

One noticeable effect is that the onset of the burst phase (ti) is delayed for higher 
metallicity and larger grain size (Table 1), as expected due to higher opacity and therefore 
slower cooling. Note that, over the bulk of our large cool disk, increasing amax increases the 
opacity. Although the time to reach the asymptotic phase is relatively insensitive to grain size 
and metallicity, the overall final tcooi hsted in Table 1 illustrates that the correlation between 
cooling time and opacity carries over to late times. During the asymptotic phase, in all cases, 
the Toomre Q values remain roughly constant with time, with values ranging between 1.3 to 
1.8 for r = 10 to 40 AU, and the mass inflow rates peak near 15 AU at ~ IQ-^Mq/yt, with 
negligible difference between 1/2 Zq and 2 Zq to the accuracy that we can measure these 
inflows (Mejia et al. 2005). Although there are some regions of superadiabatic gradients, the 
rapid overall convective cooling reported by Boss (2002, 2004) does not occur. We do see 
upward and downward motions, which we attribute to hydraulic jumps (Boley & Durisen 



- 5 - 



2006). Whether or not some of these motions are actually thermal convection, they do not 
result in rapid cooling for our disks. 

In Figure 1, which shows midplane densities at 15 ORPs, the spiral structure appears 
stronger for 1/4 than for 2 Z©. In order to quantify differences in GI strength, we compute 
integrated Fourier amplitudes (Imamura et al. 2000) 

^ _ / pmrdrdz 
^ J pordrdz ' 

where po is the axisymmetric component of the density and pm is the amplitude of the 
cos(m0) component. Although variable, after ~ 14 ORPs, the A^'s for most m's are greater 
for 1/4 Zq than for higher Z's. To measure total nonaxisymmetry, we sum the A^^s and 
average this sum over 14.5 to 15.5 ORPs. As shown in Table 1, this global measure (A) is 
greatest for 1/4 Zq and generally decreases with increasing metallicity and grain size. 

Figure 2 plots the cumulative energy loss due to coohng computed for only half the 
disk as a function of time. The upper curves show energy loss from the disk interior after 
compensating for energy input by residual irradiation and by the glowing disk upper atmo- 
sphere; the lower curves show net energy loss from optically thin regions after accounting for 
heating due to envelope irradiation. Due to our restricted vertical resolution and use of the 
Eddington atmosphere fit over one or two vertical cells, the "thick" curves effectively include 
most of the photospheric layers for most columns through the disk. The "thin" curves tally 
additional coohng from extended layers above the photospheric cells, usually with tr « 1, 
plus the parts of the outer disk that are optically thin all the way to the midplane. The 
initial coohng rates for the optically thick regions plus photosphere clearly differ. In fact, the 
initial slopes of the "thick" curves give tcooi ~ Z/Z© ORPs. However, the initial disks are far 
from radiatively relaxed, and so there are transients. Remarkably, by the asymptotic phase, 
all the disk interior-plus-photosphere curves converge to similar energy loss rates. During 
these late times, the differences between the total cooling rates are dominated by the opti- 
cally thin regions, which are larger for the lower metallicity cases, as indicated by the Thin% 
entry in Table 1. The overall asymptotic phase tcooz's in Table 1, based on summing the thick 
and thin loss rates, are longer for higher metallicity and larger grain size. Altogether, the 
evidence in Table 1 and Figures 1 and 2 shows that higher opacity leads to slower cooling 
and that slower cooling produces lower GI amplitudes. We remind the reader that we detect 
these differences over a much narrower range of metallicities (1/4 to 2 Zq) than considered 
by Boss (2002) (0.1 to 10 Zq). 

As in Mejia (2004), except for brief transients during the burst phases of some runs, these 
disks do not form dense clumps, in apparent disagreement with Boss (2002). To investigate 
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whether the disk evolution depends on spatial resolution in the asymptotic phase (Boss 2000, 
2005; Pickett et al. 2003), both the 1/4 and 2 Zq simulations are extended for another 
2 ORP's with quadrupled azimuthal resolution (512 zones), and the disks do not fragment 
into dense clumps. This is consistent with the analytic arguments in Rafikov (2005) that 
an unstable disk and fast radiative cooling are incompatible constraints for realistic disks at 
10 AU (see Boss 2005 for a different perspective). Indeed, if tcooi listed in Table 1 is a good 
measure of local coohng times in these disks, we do not expect fragmentation. Gammie 
(2001) shows that fragmentation occurs only if the local tcooi is less than about half the 
local disk orbit period Prot (see also Rice et al. 2003; Mejia et al. 2005), except possibly 
near sharp opacity edges (Johnson & Gammie 2003). We only find locallized cooling times 
shorter than 0.5 Prot in the asymptotic phase of the 2 Zq case, and then only in the 30 to 40 
AU region, which is optically thin. This occurs because, even though tcooi ~ Z in optically 
thick regions (higher optical depth), tcooi ~ in thin ones (more emitters). As a result, 
this disk displays the steepest drop of local tcooi with r. The short local tcooi' s appear to 
be highly variable and transient. The continuation of this simulation for 2 ORPs at higher 
azimuthal resolution (512 zones) does not show evidence for fragmentation into clumps. It 
could prove important to push our simulations to higher Z in the future. 



4. DISCUSSION 

Our results show that GI strength decreases as metallicity increases and, contrary to 
Boss (2002), that global radiative cooling is too slow for fragmentation into dense clumps. 
In the asymptotic phase, cooling rates for the disk interior plus photospheric layers converge 
for all Z, but the total coohng, including the optically thin regions, is higher for lower Z. 
Thus, the optically thin upper atmosphere and outer disk play a role in the cooling rate of an 
evolved disk. In fact, the fractional volume of the optically thin regions becomes very large 
at late times (see Thin% in Table 1). Note also that the optically thick region fractional 
volume, 1 — Thin%, varies roughly as Z. The greater surface area of the disk photosphere at 
higher Z tends to compensate for the higher opacity and leads to convergence of the cooling 
rates for the parts of the disk contained within the photospheric layers. In this respect, we 
confirm Boss's conclusion that the outcome of the radiative evolution is somewhat insensitive 
to metallicity. However, the important difference is that we do not see fragmentation into 
dense clumps, presumably because our cooling rates are much lower than in Boss (2002). 
For the 1mm case, the optically thin regions have a much smaller volume (Table 1) and 
contribute little to cooling. Outside the inner few AU, bigger grains make the disk more 
opaque to longer wavelengths, and tcooi is thus considerably larger, even initially. 
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Our results argue against direct formation of gas giants by disk instability in two ways 
- the global radiative cooling times seem too long for fragmentation to occur and GI's are 
stronger overall for lower metallicity. Nevertheless, it is still possible that GI's play an 
important role in gas giant planet formation. Durisen et al. (2005) suggest that dense gas 
rings produced by GI's will enhance the growth rate of solid cores by drawing solids toward 
their centers (Haghighipour & Boss 2003) and thereby accelerating core accretion. Such 
rings are indeed produced in the inner disks of all our calculations regardless of metaUicity 
or grain size, and they appear to be still growing when the calculations end. In the weaker 
GI environments of high metallicity, there is less self-gravitating turbulence to interfere with 
the radial drift of solids (Durisen et al. 2005). In this way, rings may provide a natural 
shelter and gathering place for growing embryos and cores. 

The apparent disagreement between our results and those of Boss (2002, 2004) could be 
due to any number of differences in techniques and assumptions, such as artificial viscosity, 
opacities, equations of state, initial disk models and perturbations, grid shapes and resolu- 
tion, and radiative boundary conditions, including the way that we handle irradiation. We 
are now collaborating with Boss in an effort to pinpoint which of these is the principal cause 
(K. Cai et al., in preparation). Preliminary results suggest that it is the radiative boundary 
conditions. We are therefore developing alternative techniques for disk radiative transfer 
that we hope are more reliable and accurate. 

Wc thank A. P. Boss and an anonymous referee for useful comments. This work was sup- 
ported in part by NASA Origins of Solar Systems grants Nos. NAG5-11964 and NNG05GN11G, 
by NASA Planetary Geology and Geophysics grant No. NAG5-10262, and by a Shared Uni- 
versity Research grant from IBM, Inc. to Indiana University. 
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Table 1: Simulation Results 



Case 


Ik 




Duration^ 






(A) 


tcool 


Thin% 


1/4 Ze 


1/4 


1 jim 


3.8^ 


N/A 


N/A 


1.29 


2.1 


99% 


1/2 Zo 


1/2 


1 fim 


15.6 


4.0 


10 


1.09 


2.9 


98% 


Ze 


1.0 


1 fim 


15.7 


5.0 


10 


1.10 


3.2 


94% 


2Zq 


2.0 


1 fim 


16.5 


5.0 


10 


0.72 


3.7 


86% 


1mm 


1.0 


1 mm 


17.2 


7.0 


11 


0.88 


4.5 


44% 



"AH times are given in units of ORPs. 
''Prom 13.0-16.8 (ORPs). 



Fig. 1. — Midplane density maps at 15 ORPs for the 1/4 (left panel) and 2 (right 
panel) simulations. Each square is 113 AU on a side. Densities are displayed on a logarithmic 
scale running from light grey to black (print version) or dark blue to dark red (online version) , 
as densities range from about 4.8x10"^^ to 4.8x10"" g cm , respectively, except that both 
scales saturate to white at even higher densities. 
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t (ORP) 

Fig. 2. — Cumulative total energy loss as a function of time due to radiative cooling in 
optically thick (upper set, labelled "thick") and optically thin regions (lower set, labelled 
"thin"). Both of these are net global cooling after heating by irradiation is subtracted. The 
curves labeled by a metaUicity value all use Qmax = 1 A^m. The curves labeled "1mm" are for 
a calculation with amax = 1mm and solar metallicity. Note that the 1/4 run starts from 
the 1/2 Zq simulation at about 13 ORPs. A color version of this figure appears on line. 



